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Abstract 



The ongoing discovery of terrestrial exoplanets accentuates the importance of studying planetary evolution for 
a wide range of initial conditions. We perform thermal evolution simulations for generic terrestrial planets with 
masses ranging from that of Mars to lOM® in the stagnant-lid regime, the most natural mode of convection with 
strongly temperature-dependent viscosity. Given considerable uncertainty surrounding the dependency of mantle 
rheology on pressure, we choose to focus on the end-member case of pressure-independent potential viscosity, where 
viscosity does not change with depth along an adiabatic temperature gradient. We employ principal component 
analysis and linear regression to capture the first-order systematics of possible evolutionary scenarios from a large 
number of simulation runs. With increased planetary mass, crustal thickness and the degree of mantle processing 
are both predicted to decrease, and such size effects can also be derived with simple scaling analyses. The likelihood 
of plate tectonics is quantified using a mantle rheology that takes into account both ductile and brittle deformation 
mechanisms. Confirming earlier scaling analyses, the effects of lithosphere hydration dominate the effects of planetary 
mass. The possibility of basalt-eclogite phase transition in the planetary crust is found to increase with planetary 
mass, and we suggest that massive terrestrial planets may escape the stagnant-lid regime through the formation of 
a self-destabilizing dense eclogite layer. 
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1. Introduction 

Plate tectonics is only observed on Earth and is likely 
important to Earth's uniquely clement surface condi- 
tions (e.g., Kasting and Catling[ 20031. Other terres- 
trial planets in the Solar System (i.e.. Mercury, Mars, 
and Venus) are generally considered to feature a rigid 
spherical shell encompassing the entire planet, with 



hot mantle convecting beneath the shell (e.g., Schubert 



et al. 2001 1. This mode of mantle convection is known 
as stagnant-lid convection. In fact, stagnant-lid con- 
vection may be most natural for planetary mantles be- 
cause the viscosity of constituent materials is strongly 



temperature-dependent (Solomatov, 19951. The discov- 



ery of many extrasolar terrestrial planets with mass 1 to 



lOMe (e.g.. 


Rivera et al. 


2005 


Udry et al. 


2007 


Queloz 


et al.j 20091 


Mayor et al. 


2009 


Leger et al 


2009 


Char- 


bonneau et al. 2009 Borucki et al.| 2011 ) 


makes under- 



standing planetary evolution in the stagnant-lid regime 
especially critical. 

Parametrized models of stagnant-lid convection have 
long been applied to planets in our Solar System in 
an effort to infer likely planetary evolution scenarios 



from limited observational constraints (e.g.. 


Stevenson 


et al. 


1983 


SpohnJ 


1991 


Hauck and Phil 


ipsl 2002| 


Fraeman and Korenaga 


2010 


). Previous studies of 



massive terrestrial planets are more theoretical in na- 
ture, focusing on two broad questions. First, the effects 



of planetary mass on the likelihood of plate tectonics 
have been studied through scaling analyses and simple 



parametrized convection models (Valencia et al. 2007 
O'Neill and Lenardic[ [20071 |Korenaga[|2010a[ [van Heck| 



and Tackley 2011). Second, the evolution of planets 
in the stagnant-lid regime has been contrasted with 
evolution with plate tectonics in the hope of identi- 
fying atmospheric signatures that would indicate the 
regime of mantle convection for a distant planet (e.g.. 
Kite et al. 2009). Mantle dynamics in the stagnant-lid 



regime, however, can be more complex than previously 
thought owing to the effects of mantle processing and 
crustal formation, and the scaling law of stagnant-lid 
convection that takes such complications into account 
has been developed only recently (Korenaga 2009). It 
is thus warranted to take a fresh look at the fate of mas- 
sive terrestrial planets in the stagnant-lid regime and to 
explore the general effects of initial conditions including 
planetary mass. 

This study extends a parametrized model of stagnant- 



lid convection recently applied to Mars (Fraeman and 



Korenaga, 2010) to terrestrial planets of various masses, 
including massive planets that evolve in the stagnant- 
lid regime that are termed "super- Venus" planets. This 
model incorporates the effects of compositional buoy- 
ancy and dehydration stiffening on mantle dynamics 



(Korenaga 2009j), which are rarely accounted for eyi- 
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cept in simulations of plate tectonics. Unlike in previous 
studies, sensitivity analyses are extensively performed 
to quantify the relationship between initial conditions 
and modeling results. Principal component analysis is 
used to simplify the interpretation of a large number of 
simulation results. Simple scaling analyses are also con- 
ducted to derive a theoretical basis for major modeling 
results. Moreover, the likelihood of plate tectonics is 
quantified by tracking the viscosity contrast across the 
lithosphere during each simulation. 

The purpose of this study is to investigate paths along 
which generic terrestrial planets may evolve and to esti- 
mate whether massive terrestrial planets are relatively 
more or less likely to escape the stagnant-lid regime. 
Throughout this paper, "Mars" and "Venus" should 
be considered shorthand for generalized O.IOTM® and 
O.8I5M0 terrestrial planets, respectively. The evolu- 
tion of a particular planet is likely to diverge from the 
predictions of these simple parametrized models. Few 
constraints are available beyond planetary mass and ra- 
dius for extrasolar terrestrial planets. But for terres- 
trial planets in our Solar System, more data are avail- 
able from decades of observations and spacecraft visits. 
Here, we explore hypothetical planetary evolution with 
the simplest assumptions on mantle dynamics, thereby 
serving as a reference model on which additional com- 
plications may be considered if necessary. 

2. Theoretical Formulation 

Parametrized convection models are used to simulate 
the evolution of Mars, Venus, and putative super- Venus 
planets for a wide range of initial conditions. Equa- 
tions used to track the thermal and chemical evolu- 
tion of terrestrial planets are taken from [Fraeman and| 
Korenaga (2010) with some modifications. Earth-like, 



peridotite mantle compositions are used to parametrize 
melting behavior. Although continuous evolution in the 
stagnant-lid regime is assumed, a simple model of litho- 
spheric weakening is also considered to evaluate the like- 
lihood of plate tectonics occurring at some point during 
planetary evolution. 

2.1. Governing Equations 

Mars, Venus, and super- Venus planets are assumed 
to begin as differentiated bodies with a mantle and 
core. Energy conservation yields two governing equa- 



tions. First, the energy balance for the core is ( Steven- 



son et al. 1983 1 



dRi An 



dt 



iTTRiF,, 

(1) 

where i?c and Ri are the radii of the core and inner core, 
respectively; pc is the density of the core; Lc is the latent 
heat of solidification associated with the inner core; Eg 
is the gravitational energy liberated per unit mass of the 
inner core; 77c is the ratio of Tcm, the temperature at the 
core side of the core/mantle boundary, to the average 



core temperature; Cc is the specific heat of the core; 
and Fc is the heat fiux out of the core. The formulation 
of core cooling is identical to that of [Stevenson et ah] 
p83|. 



Second, the energy balance for the mantle is (Hauck 



and Phillips 2002 ) 



An 



dt 

MRrnFn 



Prn fm^m 

-RlFc), (2) 



where Rm is the radius of the mantle; Q,„ is the vol- 
umetric heat production of the mantle; pm is the den- 
sity of the mantle; Cm is the specific heat of the man- 
tle; rjm is the ratio of the average temperature of the 
mantle to T„, the potential temperature of the man- 
tle (a hypothetical temperature of the mantle adiabati- 
cally brought up to the surface without melting); is 
volumetric melt production with associated latent heat 
release, L™; and Em is the heat flux across the man- 
tle/crust boundary. 

Some of the above parameters are universal con- 
stants, but most are planet-specific. Many important 
parameters are also time-varying. In particular, mantle 
melt is extracted to form crust, causing Rm to decrease 
with time. Likewise, Qm decreases with time because of 
radioactive decay with some approximated average de- 



cay constant, A (Stevenson et al. 19831, and extraction 



through mantle processing. 

Figure [T] illustrates the assumed thermal and chem- 
ical structure in our model. Over time, melting pro- 
cesses an upper region of the original primitive man- 
tle (PM) to form the crust and the depleted mantle 
lithosphere (DML). In parallel, the mantle lithosphere 
(ML), which is always thicker than the DML, develops 
as a conductive thermal boundary layer underlying the 
crust. As part of the DML can potentially delaminate 
and be mixed with the convecting mantle, the composi- 
tion of the convecting mantle can be more depleted than 
that of the PM. The mantle below the DML is thus re- 
ferred to as the source mantle (SM), the composition of 
which is initially identical to the composition of the PM 
but can deviate with time. The history of these layers 
strongly depends on convective vigor, effects of mantle 
melting, and initial conditions. 

2.2. Stagnant-Lid Convection with Mantle Melting 

Standard parameterizations are used for mantle rhe- 
ology and the vigor of convection. Mantle viscosity is 
a function of mantle potential temperature and the de- 



gree of hydration as (Fraeman and Korenaga 2010) 



niTu.c, 



W \ 
SM) 



A exp 



E 
RXu 



(3) 

where A is a constant factor calculated using a reference 
viscosity 770 at a reference temperature T* — 1573 K; 
E is the activation energy; R is the universal gas con- 
stant; and At/^ is the viscosity contrast between wet 
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Figure 1: Cartoons showing the assumed thermal and chemical structure of terrestrial planets taken from |Fraeman and Korenaga||2010[ |. 
In general, terrestrial planets are divided into a crust, mantle, and core, as shown in the left panel. Mantle that has been processed 
by melting and stays in the thermal boundary layer is depleted mantle lithosphere (DML). The thickness of the DML must always 
be equal to or less than the thickness of the mantle lithosphere (ML). The section of the mantle below the thermal boundary layer is 
the sublithospheric mantle. The right panel shows the horizontally- averaged temperature distribution. Key model parameters are also 
indicated. 



and dry mantle. We use an activation energy of E = 
300 kJ mol^^, which is appropriate for diffusion creep 
and dislocation creep within a Newtonian approxima- 



naga 



tion (Christensen 1984 Karato and Wu 1993 Kore- 



2006). The normalized water concentration in 
the source mantle, C^^, has an initial value of one 
and decreases toward zero as mantle melting causes 
dehydration. To write Eq. |3j we make the major as- 
sumption that mantle viscosity is not strongly pressure- 
dependent, which is consistent with early studies of the 



evolution of large rocky planets (e.g., Valencia et al. 



2006) and some theoretical predictions (Karato 2011) 



but in contrast to recent work (e.g., Papuc and Davies 



2008[ [Stamenkovic et al.[ |2011[ |20i2[ ). Our choice is 



thus further explained in the discussion section. 

Two non-dimensional parameters characterize ther- 
mal convection with the above viscosity formulation 



(Solomatov 1995). First, the internal Rayleigh number 



serves to quantify potential convective vigor (Fraeman 



and Korenaga 2010 ) 



i?a, — 



K?y(T„, C^) 



(4) 



where a is the coefficient of thermal expansion; k is the 
thermal diffusivity; Tc and are, respectively, the tem- 
perature at the bottom of the crust (called "Moho tem- 



perature") and the mantle potential temperature de- 
fined at the top of the mantle; and hm is the thickness of 
the mantle. Second, the Frank-Kamenetskii parameter 
is defined as ( Solomatov] 1995[ Fraeman and Korenaga 
20T0| 



RT?. 



(5) 



With these two parameters, the average convective 
velocity beneath the stagnant-lid may be calculated as 



( Solomatov and Moresi 2000 ) 



0.38- 



Ra 



1/2 



(6) 



To include the effects of compositional buoyancy and 
dehydration stiffening, the Nusselt number, which is a 
non-dimensional measure of convective heat flux, must 
be calculated with a local stability analysis at each time 



step ( Korenaga 2009 ) . The symbolic functionality may 



be expressed as 



Nu = f{Ra,E,Tu,Tc,hi,h^,Ap,Ar]„ 



(7) 



where Ap and Arjm are the density and viscosity con- 
trasts between the source mantle and depleted mantle, 
respectively, and hi is the thickness of the depleted man- 



3 



tie lithosphere. The thickness of a thermal boundary 
layer in the mantle is then easily calculated using 



'ML 



(8) 



The chemical evolution of the mantle strongly affects 
terrestrial planet evolution. To first order, partial melt- 
ing of the mantle can be considered to begin at a depth 
where the temperature exceeds the solidus of dry peri- 
dotite, as long as the mantle is not significantly wet 
(Hirth and Kohlstedt 1996). The initial pressure of 



melting is (Korenaga et al. 2002) 



Pi = 



Tu- 
1.20 X 10-7 



1423 

-idT/dP)s' 



(9) 



where {dT/dP)s is the adiabatic mantle gradient, which 
is roughly constant for the pressure range relevant to 
mantle melting. Therefore, Pi should be approximately 
constant for any terrestrial planet with Earth-like man- 
tle composition. Melting stops when the convective up- 
welling reaches the base of the mantle lithosphere. That 
is, the final pressure of melting is given by 



Pf = PLg{K + Jt-ml) 



(10) 



where is the thicknesses of the crust; g is gravita- 
tional acceleration; and is the density of the litho- 
sphere. For convenience, we use the Martian pm as pi, 
noting that should remain roughly constant whereas 
Pm, an averaged mantle parameter, increases with plan- 
etary mass because of pressure effects. If Pf < Pi, then 
melting occurs in the melting zone between Pi and Pf, 
with thickness dm and average melt fraction equal to 



P-P. 



d<p 
dP 



(11) 



where {d(f>/dP)s is the melt productivity by adiabatic 
decompression. Volumetric melt productivity is finally 
parametrized as 



I. 



2xdrnU(j) 2 
; 47rit^ 



(12) 



where x ~ 1 if the upwelling mantle is cylindrical and 



all downwelling occurs at the cylinder's edge ( 


Soloma- 


tov and Moresi 


2000 


Fraeman and Korenaga 2010). 


The crustal temperature profile is calculated as in Frae-| 



man and Korenaga (2010), with the modification that 
crustal material with a temperature above Tcr it = 1273 
K is considered to be buoyant melt that migrates within 
one time step immediately below the planet's surface, 
producing a relatively cooler crust and a larger man- 
tle heat flux. This modification is not important for 
Martian cases, because the Moho temperature does not 
reach the threshold except for some extreme cases, but 
becomes essential to achieve a realistic crustal thermal 
profile for larger planets. 



2. 3. Likelihood of Plate Tectonics 

Thermal evolution models featuring stagnant-lid con- 
vection are not applicable to planets on which plate 
tectonics occurs. If a suitable weakening mechanism 
exists, the lithosphere may be broken into plates and 
recycled into the mantle. Many aspects of plate tec- 
tonics on Earth, however, are not captured in current 
mathematical models (Bercovici 2003). Quantifying 



the conditions under which plate tectonics is favored 
over stagnant-lid convection is likewise difficult, and the 
effect of planetary mass on the likelihood of plate tec- 



tonics has been controversial ( 


Valencia et al. |2007 Ko- 


renaga 




2010a 




van Heck and Tackley 




2011|). Recent 



studies suggest, however, that the effects of planetary 
mass on yield and convective stresses may be dominated 
by uncertainties in other important planetary parame- 
ters, such as internal heating and lithosphere hydration 



(Korenaga 2010a van Heck and Tackley 2011) 



This study uses a simple scaling that is consistent 



with current understanding of rock mechanics (Kore- 



naga 2010a), though the possibility of different litho- 



sphere weakening mechanisms (e.g., Landuyt et al 



2008) cannot be excluded. We assume that plate tec- 
tonics can occur if convective stress exceeds the brittle 
strength of lithosphere given by 

Ty = co + ppgz, (13) 

where cq is the cohesive strength, p is the effective fric- 



tion coefficient, and z is depth (Moresi and Solomatov 



1998). Experimental data indicate that the cohesive 



strength is negligible under lithospheric conditions, i.e. 



co/(/xpz) ^ 1 (Byerlee 19781. We use another non- 
dimensional parameter (Korenaga 2010a I: 



7 



(14) 



where the relevant temperature difference is the differ- 
ence between the mantle potential and surface temper- 
atures. In the parameterized convection model formu- 
lated in the previous section, we separately consider the 
crust and the mantle, but when discussing the likeli- 
hood of plate tectonics using the scaling of |Korenaga| 
(2010b), it is more convenient to treat the crust and 



mantle together, assuming that crustal rheology is sim- 
ilar to mantle rheology. 

Detailed scaling analyses (Korenaga 2010a|b ) show 
that the effective viscosity contrast across the litho- 
sphere can be parameterized as 



A?7i = exp(O.3277"''^'0tot), 



(15) 



where 6tot is the Frank-Kamenetskii parameter defined 
using the total temperature difference explained above, 
i.e., 

EiTu - T.) 



i?T2 



(16) 



A transition from plate-tectonic to stagnant-lid convec- 
tion can take place if the above viscosity contrast ex- 



4 



ceeds a critical value 



A77 



L.crit 



= 0.25i?a 



1/2 

i,tot^ 



(17) 



where Rai^tot is defined to incorporate surface temper- 
ature as 



Ra, 



apmg{Tu - Ts){hc + hmf 



(18) 



For each simulation, if Ai]]^/ Arji^ crit ^ 1 at any time, 
then plate tectonics may have been favored at some 
point during the evolution of a given planet. The sat- 
isfaction of this criterion may strongly depend on the 
value of /z, so a wide range of values should be tested. 
For silicate rocks, plausible values of /x range from 0.6 



to 0.7, according to both laboratory studies (Byerlee 



1978) and measurements of crustal strength on Earth 



e.g., Brudy et al. 19971. Surface water, however, may 



lower these values substantially via thermal cracking 



and mantle hydration (Korenaga 2007). 



Constant 


Value 


T T * 4-^ 

Units 


Ket. 


A 


1.38 X 10^^'' 


s ^ 


[1] 


k 


4.0 


W m"^ 


[1] 


a 


2 X 10"*' 


K ^ 


[1] 


K 


lO"*" 


9 1 
m s 


[1] 


PL 


3527 


1 —3 

kg m ^ 


ri 1 


Lm 


6.0 X 10^ 


Jkg-i 


[2] 


Lc + Eg 


1.0 X 10~6 


T 1 —1 

J kg 


[2] 




1 Anna 
iOOU 


T K.^ — 1 — 1 

J Kg K 




Cc 


850"^ 


J kg~i 


[3] 


Vm 


1.3^^ 


N/A 


[1] 


Vc 


1.2^^ 


N/A 
K Pa-i 


[1] 


{dT/dP)s 


1.54 X 10-* 


[4] 


{d^/dP)s 


1.20 X 10-* 


Pa~i 


[4] 



Table 1: Summary of universal constants used in all simulations. 
Reference s: 1. [Stevenson et al.| | |1983[l, 2. [Fraeman and Korenaga] 
l |2010[ ), 3. |Noack et al.| l |2011[ |, 4. |Korena!ga et al.| | |2002||. "Mars 
has Crn = 1149; Cc = 571; Tjm = 1-0; and rjc = 1.1 l |Fraeman and| 
|Korenaga||2010| . 



3. Numerical Models 

The parametrized model described above was used 
to calculate thermal histories of Mars- and Venus-like 
planets and 1 to lOM^ super- Venus planets, where the 
subscript denotes parameters for Earth, for a du- 
ration of 4.5 Gyr using numerical integration with a 
time step of 1 Myr. A wide parameter space was ex- 
plored by varying the initial mantle potential tempera- 
ture, T„(0); the initial core/mantle boundary tempera- 
ture, T'cm(O); the initial volumetric heat production, Qg; 
the reference mantle viscosity, rjo; and the viscosity con- 
trast between dry and wet mantle, Arj^^,. Previous work 
for Mars demonstrated that simulation results were not 
very sensitive to the degree of compositional buoyancy 



inner core growth was disallowed and total surface heat 
flux at the present was required to be positive. Further- 
more, the condition hc{tp) < 500 km was imposed to 
disregard results with unrealistic crustal growth. None 
of the 675 simulations for Venus failed these criteria, 
but 30 of the 675 simulations for Mars were discarded. 

The appropriate magnitude of radiogenic heating is 
poorly constrained in general, especially for terrestrial 
exoplanets. Even for Earth, the abundance of radio- 
genic heating is controversial. Geochemical constraints 
support a low Urey ratio, the ratio of internal heat 
production to surface heat flux, but this is known to 
conflict with the cooling history of Earth unless a non- 
classical heat-flow scaling for mantle convection is as- 
sumed (Christensen 1985 Korenaga 2008). A Urey 



and other parameters ( [Fraeman and Korenaga] |2010D . ^.^tio close to one has thus long been preferred from a 



Table[T]lists model constants common to all simulations, 
and Table [2] lists planet-specific ones. 

3.1. Application to Mars- and Venus-like Planets 

For Venus, the following sets of initial conditions were 
used: Initial mantle potential temperature, Tu{0) = 
1400, 1550, 1700, 1850, and 2000 K; initial core/mantle 
boundary temperature, Tcm{0) ~ 3500, 4000, and 
4500 K; reference viscosity, 770 = 10^*, 10^^, and lO^o Pa 
s; and dehydration stiffening, A?/^ = 1, 10, and 100. 
For Mars, initial core/mantle boundary temperatures 
were 2250, 2500, and 3000 K. In addition, five different 
values were tested for the amount of internal heating 
Qq. Compositional buoyancy was set as (dp/dcf)) = 120 
kg m-'^ for all simulations. The fraction of light el- 
ements in the core was fixed at 0.2 for all simulation 



geophysical perspective (Davies 1980 Schubert et al. 



1980 2001 ) and can be used to provide an upper bound 
for the initial concentration of radioactive elements in 
Earth's chemically undifferentiated mantle. Assuming 



a present-day surface heat flux of 46 TW ( Jaupart et al 



2007), an extreme upper bound for Earth is Qo « 
3.5x10"^ W m-3. 

A recent petrological estimate on the thermal history 
of Earth is actually shown to favor a low Urey ratio 



■^0.3) with a non-classical heat-flow scaling (Herzberg 



et al. 2010), indicating that geochemical constraints on 
the heat budget may be robust. In the thermal evo- 
lution models of Kite et al. (2009), for example, con- 



centrations of 40k, ^''^Th, ^35u77nd ^^^V taken from 



Ringwood 


( 


1991 


) and 



runs to avoid inner core solidification (Schubert et al 



1992 Fraeman and Korenaga 2010). Simulations were 



performed for all permutations of the above initial con- 
ditions, although unrealistic simulation results were dis- 
carded in the following way: For both Venus and Mars, 



were considered, corresponding to values for Qq between 
1.2 X 10"'^ and 8.7 x 10"* W m'^ for Venus. We thus 
chose to use the following values for initial volumetric 
radiogenic heating: Qo — 0.5, 0.75, 1.0, 1.25, and 1.75 
xlO"*" W m~'^ (in the case of Venus). The default in- 
termediate value is 1.0 xlO~^ W m"-^. For other plan- 
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Constant 


Mars 


Venus 


lA/e 


2Me 


4Me 


5Me 


6Me 


8Me 


lOM® 


Units 


9 


3.70 


8.87 


10.0 


13.6 


18.6 


20.7 


22.6 


26.0 


29.1 


ni s 


Ts 


220 


730 


300 


300 


300 


300 


300 


300 


300 


K 


tip 


ooyO 


CCA 


OAT 


(bm 


92oz 


AO0 1 


1 AO A C 

10295 


1 1 ATO 


1 1 add 
11696 


km 


Rc 


1550 


3110 


3295 


3964 


4723 


4986 


5206 


5564 


5848 


km 


Pm 


3527 


3551 


4476 


4951 


5589 


5845 


6078 


6497 


6873 


kg m^"^ 


Pc 


7200 


12500 


12961 


14882 


17594 


18698 


19708 


21530 


23174 


kg m^"^ 


P 

^ cm 


19 


130 


151 


284 


556 


697 


842 


1144 


1463 


GPa 


Pc 


40 


290 


428 


821 


1639 


2067 


2508 


3431 


4406 


GPa 



Table 2: Summary of planet-specific constants for Mars, Venus, and seven super- Venus planets. Martian values were taken from |Fraeman| 
[and Kore naga (2010) and references therein. Venusian values can b e found in|Spolm (1991j and |Noack et al^ ( |2011^ . Super- Venus values 
were calculated in this study from simple interior models following [Seager et al.| l |2067^ 



ets, Qq was multiplied by Pm/Pm $i where $ indicates 
the Venusian value, to maintain constant element abun- 
dances in more or less compressed mantles. 

3.2. Application to Super- Venus Planets 

One-dimensional profiles of massive terrestrial exo- 
planets were generated to calculate planet-specific con- 
stants used in the above stagnant-lid convection model. 
Many interior structure models exist for massive solid 



exoplanets, ranging from simple to very complex (Va- 



lencia et al.t,2006;.Scager et al .[[200711 Sot in et al.[[2007 



Wagner et al. 2011 1 . To study the first-order effects of 



planetary mass on stagnant-lid convection, a relatively 
simple structure with an Fe(e) core and a MgSiOs man- 
tle is assumed, as in Seager et al. ( 2007 ) and Kite et al. 



(2009). The resulting interior density and pressure dis- 



tributions neglect several obvious factors such as tem- 
perature effects, but yield results remarkably similar to 
those from more complex models. 

Three equations are solved to calculate m(r), the 
mass contained within radius r; P{r), the pressure dis- 
tribution; and p{r), the density distribution. A self- 
consistent internal structure must satisfy the material 
specific equations of state 



P{r) = fEOs{p{r),T{r)), 
the equation of hydrostatic equilibrium 

dP{r) _ -Gm{r)p{r) 
dr ' 



(19) 



(20) 



and the conservation of mass equation for a spherical 
mass distribution 



dm{r) 
dr 



47rr^p(r), 



(21) 



where G is the gravitational constant; T(r) is the radial 
temperature profile; and fsos represents a material- 



specific equation of state (Seager et al.[ 2007). 



Equations of state are numerically calculated using 
constants from Table [H to sufficient resolution so that 
p{r) can be determined to within ±1 kg m^'^. For a 
desired Mp, equations (20) and (21) are numerically 




0.5 1 
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Figure 2: Interior density distributions for super- Venus planets 
with Mp = 1, 2, 4, 5, 6, 8, and lOM^g. The material equations of 
state and the equations of conservation of mass and hydrostatic 
equilibrium were numerically integrated to build simple planets, 
and the interior boundary condition was adjusted until the result- 
ing planet had the desired mass. These simple models are used to 
calculate averaged values for mantle density, p^; core density, p^; 
surface gravity, g; central pressure, Pc; and core/mantle bound- 
ary pressure, Pcm- For comparison, the density distribution for 
Earth from the Preliminary Reference Earth Model (PREM) of 
[Dziewonski and Anderson| | |1981[ l is also plotted. 



boundary conditions M{0) = and P{0) = Pc, where 
Pc is a guessed central pressure. The outer boundary 
condition is simply P{Rp) = 0. Errors associated with 
ignoring temperature effects are limited to a few per- 
cent (Seager et al. 2007). With this method, the choice 



of Pc determines the Rp at which the outer boundary 
condition is satisfied. These calculations are iterated 
with the bisection method until Pc is found such that 
m{Rp) = Mp to within 0.1%. The equation of state 
for Fe(e) is used until m(r) = 0.325Mp, mandating a 
32.5% core mass fraction. The MgSiOa perovskite to en- 
statite phase transition is assumed to occur at 23 GPa, 
although the transition pressure increases to ~25 GPa 
at -800 K (e.g., [Ita and Stixrude[ [1992|). Neglecting 



integrated from the center of a planet with the inner 



this phase transition would produce unrealistically high 
near-surface densities. 

Pressure, mass, and density distributions were cal- 
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Material 


ifo (GPa) 


K 


K (GPa-i) 


Pa (kg m-3) 


EOS 


Fe(e) 


156.2 


6.08 


N/A 


8300 


V 


MgSi03(pv) 


247 


3.97 


-0.016 


4100 


BME4 


MgSiOaM 


125 


5 


N/A 


3220 


BME3 



Table 3: Material constants used to generate interior structure models, taken from |Seager et aL] l |2007l l. Using three different equations 
of state, P{p) is calculated to high resolution for each material. The Vinet and 3rd and 4th order Birch-Murnaghan equations of state 
are abbreviated V, BME3, and BME4, respectively. 



culated for planets with Mp — 1, 2, 4, 5, 6, 8, and 
IOMq. From these models, averaged densities for the 
core and mantle were calculated, along with Pc and Pcm- 
Surface gravitational accelerations were calculated us- 
ing g = GMp/R^p for each planet. These constants 
are reported in Table [2] The density distributions for 
these planets are shown in Fig. [2j along with the density 
profile for Earth from the Preliminary Reference Earth 



Model (PREM) of |Dziewonski and Anderson] pMT ). 
Compared to PREM, the scheme for calculating internal 
structures used in this study overestimates the density 
of the core and underestimates the radius of the core 
of an Earth-mass planet, but returns Rp{M(^) « 
despite ignoring details of mineral composition, phase 
transitions, and temperature effects. In more massive 
planets, the enstatite to perovskite phase transition oc- 
curs at much shallower depths because a higher surface 
gravity causes a greater increase in pressure with depth. 
Furthermore, Pc increases much more rapidly than Pcm 
with increasing planetary mass. 

Most of the Martian and Venusian initial conditions 
can be used for super- Venus thermal evolution mod- 
els, but some must be modified appropriately. For ex- 
ample, the core/mantle temperature for super-Venus 
planets should increase along the mantle adiabatic tem- 
perature gradient: for 5 and lOM^ super- Venus plan- 
ets, initial core/mantle boundary temperatures are in- 
creased by roughly 350 and 900 K, respectively, from the 
initial conditions for Venus. These temperatures still 
correspond to a so-called "hot start," which is likely 
for terrestrial planets because of the large magnitude 
of gravitational potential energy released during accre- 
tion (Stevenson et al. 1983). Only three simulations 



for the 5M0 planet, and no simulations for the lOM® 
planet, failed the requirements on crustal thickness, sur- 
face heat flux, and inner core growth. 



4. Results 

Thermal evolution simulations were performed for 
Mars- and Venus-like planets and super- Venus planets. 
The following sections summarize the results, beginning 
with a few representative examples for Mars and Venus. 
Principal component analysis, as described in the ap- 
pendix, was applied using all simulation results to iden- 
tify major model behaviors. We also tried to quantify 
relations between input and output parameters and, de- 
spite the complexity of our model formulation, a linear 



function of initial conditions is found to reasonably ap- 
proximate many output parameters of interest. 

4-1. Sample Thermal Histories for Mars- and Venus- 
like Planets 

Sample thermal histories for Mars and Venus are 
shown in Fig. |3] These models span the entire range 
of initial radiogenic heating values with all other ini- 
tial conditions set to intermediate values. In particular, 
r„(0) = 1700 K, ^0 = 10^^ Pa s, and A/^„ = 100. For 
Venus and Mars, respectively, Tf,m(0) — 4000 and 2500 
K and Ts ~ 220 and 730 K. Initially very hot cores 
are assumed here because core segregation is expected 
to release a large amount of gravitational potential en- 
ergy. This excess heat is released into the mantle dur- 
ing the first hundred million years of planetary evolu- 
tion. Thereafter, mantle dynamics controls core cool- 
ing. Whereas internal heating has a great effect on sur- 
face heat flux, mantle temperatures only differ to within 
±200 K for the sampled range of Qq. Mars evolves with 
a consistently lower potential temperature than Venus. 
Because Mars also has a relatively shallow mantle, the 
Martian core is cooled down more efficiently. 

Crustal thickness is an important, potentially observ- 
able constraint for planetary evolution models. Mars 
and Venus, with different magnitudes of radiogenic 
heating, have very different crustal formation histories. 
Both planets start with no initial crust, but quickly pro- 
duce some through mantle melting. For Venus, Moho 
temperatures quickly reach the melting point of basalt 
for all initial internal heating choices. Crustal pro- 
duction occurs for the first ^1 Gyr of evolution, with 
thicker crust for higher internal heating. For Mars, 
crustal production is gradual and crustal temperatures 
are much lower, with increased internal heating causing 
an longer period of crustal formation and increased to- 
tal crustal production. Both Mars and Venus undergo 
substantial mantle processing, indicating that the deep 
interior serves as a significant source of endogenous wa- 
ter, especially during the first ^1.5 Gyr of their evolu- 
tion. 

4-2. Sensitivity Analyses for the Evolution of Mars- and 
Venus-like Planets 

Figures |4] and |5] summarize the results of 1320 simula- 
tions for Venus and Mars, respectively. Present-day val- 
ues for selected output parameters are plotted against 
crustal thickness for both planets. Several correlations 
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Figure 3: Sample liistories for Venus (top) and IVIars (bottom). From left to right, red curves signify Moho temperature, surface heat flow, 
crustal thickness, and normalized mantle water content. Likewise, blue curves represent mantle potential temperature, mantle heat flux, 
depleted mantle lithosphere thickness, and fraction of processed source mantle; green curves show core/mantle boundary temperature, 
core heat flux, and lithosphere thickness. Solid, dashed, and dotted lines indicate Venusian Qo = 0.5, 1.0, and 1.75 xlO~^ W m"'', 
respectively. Default initial conditions are r„(0) = 1700 K, 7?o = lO^'' Pa ■ s, Ar;„ = 100, and {dp/d(p) = 120 kg m'^. Venus and Mars 
have Tcm(O) = 4000 and 2500 K, respectively. Because crustal melting causes highly discontinuous surface heat flux, a moving average 
with a 75 Myr span was used for plotting purposes. 



are readily apparent. For Venus, thicker crust is as- 
sociated with higher Moho temperature, more mantle 
processing, higher mantle heat flux, and quicker crustal 
formation. More specifically, Moho temperature in- 
creases with crustal thickness in a linear fashion until 
he ~ 75 km, after which Moho temperatures remain 
near the critical value for basalt melting. Simulations 
with crustal melting have highly discontinuous surface 
and mantle heat fluxes, but such discontinuous nature is 
merely an artifact owing to our particular numerical im- 
plementation, so average values of Fg and Fm over the 
final 100 Myr of planetary evolution are used for all sub- 
sequent analyses. In contrast, Moho temperatures for 
Mars only approach the critical value for basalt melting 
in simulations with the thickest crust. Unlike for Venus, 
a decrease in present-day mantle heat flux accompanies 
an increase in crustal thickness for Mars. 

Principal component analysis facilitates the interpre- 
tation of the correlations between output parameters. 
For Venus, two principal components account for most 
(>65%) of the variance of the planetary parameters af- 
ter 4.5 Gyr of thermal evolution. Calculations of the 
principle components returns coefficients with values 
between -1 and 1 that are associated with each model 
output parameter. Appendix A contains a table of prin- 
cipal component basis vectors for Venus. Comparing 
numerical values of select coefficients may reveal corre- 
lations with physical explanations. Arrows representing 
the eigenvectors associated with these principal compo- 
nents are plotted in Figs. |4] and [S] These arrows indicate 
the axes along which the majority of the variance in 
the model output primarily lies. No preferred polarity 
exists for the principal component eigenvectors; plot- 
ting these arrows with a 180° rotation would be equally 
valid. 

The first principal component represents the most 



dominant correlations among present-day planetary pa- 
rameters, which are characterized mainly by the thick- 
nesses of the crust and mantle lithosphere layers, as 
they are associated with large coefficients: he (0.34), hi 
(-0.36), and /iml (-0-39). Because the sign of the coef- 
ficient for he is opposite to the sign of the other two co- 
efficients, the thicknesses of the crust and mantle litho- 
sphere are anti-correlated. In other words, thick crust is 
associated with thin depleted mantle lithosphere and a 
thin thermal boundary layer and vice versa, since prin- 
cipal components have no preferred polarity. An ini- 
tially hotter mantle produces thicker crust and thicker 
depleted lithosphere. Because mantle viscosity is lower 
for hotter mantle, however, the depleted lithosphere 
is more likely to be destabilized, resulting in a thin- 
ner lithosphere (and thus thermal boundary layer) for 
thicker crust. Other coefficients in the first principal 
component indicate the effects of crustal thickness on 
other model parameters, including the first-order cor- 
relations observed during inspection of Fig. |4] For in- 
stance, thick crust is associated with high Moho tem- 
perature and high surface and mantle heat fluxes. Thick 
crust also indicates a high degree of mantle processing 
and a corresponding low present-day mantle water con- 
tent. Finally, the large negative coefficients for both 
tc,io% s-'^d log]^Q(Aic, tot) indicate that thick crust tends 
to form early and quickly. 

The second principal component elucidates the ef- 
fects of planet temperatures on other model parameters, 
since large coefficients are associated with (0.45) and 
Tcm (0.44). Unsurprisingly, high mantle potential and 
core/mantle boundary temperatures are associated with 
high Moho temperature, since Tc has a coefficient of 
0.25. Moreover, high interior temperatures correspond 
to thick crust and a high degree of mantle processing, 
which would cause the present-day mantle water con- 
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Figure 4: Summary of parameter values at the present for 675 simulations of the thermal evolution of Venus. Arrows are projections 
of the principal component basis vectors that emanate from a point representing the averaged simulation results, indicating axes that 
account for the vast majority of the data set's variance. The red arrow represents a larger percentage of cumulative variance (41%) 
than the green arrow (27%). Panels show (a) Moho temperature, (b) mantle potential temperature, (c) surface heat flux, (d) mantle 
heat flux, (e) fraction of mantle processed by melting, and (f) total time for crust to grow from 10% to 95% of its present thickness as 
functions of crustal thickness. Because crustal melting causes highly discontinuous surface and mantle heat fluxes, the model outputs 
are the averaged values for the final 100 Myr of planetary evolution. 



centration to be very low. Note, however, that high 
(present-day) interior temperatures do not correspond 
to thick crust in case of the first principal component 
(Fig. 4b). In this space of crustal thickness and upper 
mantle temperature, the first and second principal com- 
ponents are nearly orthogonal, thus explaining the over- 
all spread of simulation results. With principal compo- 
nent analysis, we can visualize how the most dominant 
trend (represented by the first principal component) is 
affected by secondary factors and how these secondary 
factors manifest in different parameter spaces. An im- 
portant point is that the overall variability of plane- 
tary evolution can be compactly represented by a small 
number of principal components; that is, the effective 
dimension of the model space is actually small. 

The principal components for Mars are very similar 
to those for Venus, with some notable exceptions. The 
first principal component again represents the effects 
of strongly correlated Moho temperature and crustal 
thickness, and thus explains the largest portion of the 
variance in the model output. As for Venus, a thin, 
cold crust is associated with thick depleted mantle litho- 
sphere, a thick thermal boundary layer, a low surface 
heat flux, and a low degree of mantle processing. Un- 
like Venus, however, the surface and mantle heat fluxes 
in the first principal component are anti-correlated (see 
also Fig. [5]). 

Despite the complexity of our thermal evolution 
model, some present-day parameters are found to be 



predicted with reasonable accuracy for Venus and Mars 
using a linear function. A general formula for this func- 
tion is 

^»,3(logio(?7o))n + A:,4(logio(A?7^))„ + 5^0,71, (22) 

where Bi is the value of the desired output parame- 
ter after 4.5 Gyr, constants A^^o through Ai^^ are esti- 
mated using the least-squares method for each Bi, and 
the subscript n indicates that the input parameters are 
normalized and mean subtracted. 

Table |4] lists constants for Venusian Bi that have rel- 
atively high correlation coefficients between predicted 
and actual simulation results. Figure [6] shows contour 
plots with predicted values of mantle potential temper- 
ature, crustal thickness, and duration of crustal forma- 
tion for given initial internal heating and mantle poten- 
tial temperature. While present-day mantle potential 
temperature and crustal thickness depend strongly on 
both initial mantle potential temperature and the mag- 
nitude of internal heating, the total duration of crustal 
formation is primarily a function of initial mantle poten- 
tial temperature (see Table |4] for more complete infor- 
mation on parameter sensitivity). Figure |6] also demon- 
strates a reasonable correspondence between the pre- 
dicted and actual values of these model outputs for all 
of the simulations. This way of summarizing simula- 
tion results allows us not only to see the sensitivity of 
model outputs to initial parameters but also to quickly 
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Figure 5: Summary of parameter values at the present for 645 simulations of the thermal evolution of Mars. Arrows are projections 
of the principal component basis vectors that emanate from a point representing the averaged simulation results, indicating axes that 
account for the vast majority of the data set's variance. The red arrow represents a larger percentage of cumulative variance (42%) 
than the green arrow (23%). Panels show (a) Moho temperature, (b) mantle potential temperature, (c) surface heat flux, (d) mantle 
heat flux, (e) fraction of mantle processed by melting, and (f) total time for crust to grow from 10% to 95% of its present thickness as 
functions of crustal thickness. Because crustal melting causes highly discontinuous surface and mantle heat fluxes, the model outputs 
are the averaged values for the final 100 Myr of planetary evolution. 



reproduce major modeling results without redoing sim- 
ulation. 

4-3. Evolution of Super- Venus Planets 

We investigate the evolution of super- Venus planets 
to explore the effects of planetary mass on stagnant-lid 
convection. For simplicity, surface temperatures for all 
super- Venus planets are assumed to be 300 K, though 
in reality this temperature may vary with time and is 
highly dependent on atmospheric composition and on 
the distance to and the luminosity of the central star. 

4.3.1. Sample Thermal Histories 

Super- Venus planets with Mp = 1,5, and lOM® were 
evolved to study the effects of increasing planetary mass 
on a variety of parameters, particularly crustal produc- 
tion. For all three planets, Qq was scaled to the Venu- 
sian value of l.OxlO"'^ W m'^, T„(0) = 1700 K, and 
r„„(0) = 4000, 4350, and 4900 K, respectively. De- 
hydration stiffening and compositional buoyancy were 
both incorporated as usual. Figure [7] shows the results 
of these simulations. As with Venus and Mars, the tran- 
sient "hot start" in the core is lost in the first ^100 Myr. 
After this initial cooling, mantle dynamics controls core 
cooling. Because the mantle heats up for the first ^1 
Gyr and then cools only very slowly, core cooling is 
precluded for the first ~2 Gyr. As suggested by simple 



scaling laws (Stevenson 2003), mantle cooling paths for 



massive super- Venus planets are roughly parallel. 



Figure [7] also shows how the thicknesses of the crust, 
mantle lithosphere, and depleted mantle lithosphere 
vary with time. With increasing planetary mass, crustal 
thickness decreases. The simple scaling analyses below 
indicate that more massive planets have greater melt 
production. The observed increase in mantle potential 
temperature with planetary mass only accentuates this 
effect. The increased melt volume, however, is not suf- 
ficient to create a thicker crust on a larger planet. The 
lAi® planet in the stagnant-lid regime ceases crustal 
production soon after 1 Gyr as mantle potential tem- 
perature drops below a critical value. The increased 
interior temperatures for the more massive planets al- 
low longer durations of crustal production. For the first 
^2 Gyr of thermal evolution, the thickness of the de- 
pleted mantle lithosphere is close to that of the mantle 
thermal lithosphere, reflecting the continuous delamina- 
tion of excess depleted mantle lithosphere. Decreased 
crustal production with increasing planetary mass cor- 
responds to a smaller degree of mantle processing and 
a higher content of residual mantle water. 

4.3.2. Sensitivity Analyses 

The output of 1347 simulations for 5 and lOM^ 
super- Venus planets are shown in Figs. |8] and |9] Three 
simulations for the super- Venus planet were ex- 

cluded because they did not meet the requirements that 
he < 500 km and that inner core growth did not oc- 
cur. As for Mars and Venus, present-day parameters of 
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Figure 6: Top panels show predicted values for Venus after 4.5 Gyr given initial volumetric radiogenic heating and initial mantle potential 
temperature, the initial conditions to which the output is most sensitive. Bottom panels show the correspondence between predicted and 
actual simulation results. The dashed line represents perfect predictive power. Default initial conditions are Tcm.{0) = 4000 K, rjo = 10^^ 
Pa - s, Arjw = 100, and (dp/d(f>) = 120 kg m~'^. Panels show predicted values of (a) mantle potential temperature, (b) crustal thickness, 
and (c) total time for crust to grow from 10% to 95% of its present thickness. 



interest are plotted against present-day crustal thick- 
ness. The principal component eigenvectors, explained 
below, are projected onto each plot, emanating from the 
average simulation output. The table in the appendix 
contains the principal component basis vectors for the 
lOM® planet. 

These scatter plots reveal similarities between the 
evolution of both massive planets. For instance, Moho 
temperature increases with crustal thickness in a linear 
fashion before reaching the critical value for basalt melt- 
ing. With increasing planetary mass, the critical crustal 
thickness at which this transition occurs decreases. For 
relatively thick crust, Moho temperatures remain near 
the critical value for basalt melting. For both super- 
Venus planets, an increase in crustal thickness is asso- 
ciated with an increase in present-day mantle potential 
temperature, mantle heat flux, and degree of mantle 
processing. The total duration of crustal formation de- 
creases with increasing present-day crustal thickness. 
Again, correlations between model parameters may be 
studied in more detail with principal component analy- 
sis. 

For both planets, as for Venus and Mars, the first 
principal component is characterized by a strong cor- 
relation between Moho temperature and crustal thick- 
ness, explaining the general trends observed in Figs. |8] 
and[9j A decrease in both quantities is associated with 
an increase in the thicknesses of the depleted mantle 
lithosphere and the thermal boundary layer, a decrease 
in surface and mantle heat fluxes, an increase in the 



duration of crustal formation, and a decrease in the 
degree of mantle processing. The second principal com- 
ponent illuminates the effect of correlated interior tem- 
peratures. As expected, increasing mantle potential and 
core/mantle boundary temperatures causes an increases 
in crustal thickness, the total duration of crustal forma- 
tion, and the degree of mantle processing. 

Many present-day model parameters of interest can 
be represented as a linear function of initial conditions 
(Table [5]). Compared to the case of Venus, a greater 
number of parameters are found to be approximated 
reasonably well by this approach. The effects of melt- 
ing at the base of the crust undoubtedly remain a large 
source of nonlinearity in the model output for all terres- 
trial planets more massive than Mars. A more elabo- 
rate numerical implementation to deal with exceedingly 
high crustal temperatures may reduce such nonlinearity, 
though we did not explore this possibility. 

4-4- Scaling of Crustal Thickness and Mantle Process- 
ing 

We conduct simple scaling analyses to better under- 
stand the cause of decreasing crustal thickness and a 
decreasing degree of mantle processing with increasing 
planetary mass. 

4.4- 1- Crustal Thickness 

A number of parameters govern the scaling of crustal 
thickness with planetary mass. Increased melt produc- 
tion, for instance, is the flrst requirement for thicker 
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K 
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he 
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26.4 


22.2 
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-7.84 
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0.92 


hi 


53.6 


-2.17 
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Fs 
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p 


29.6 


3.40 
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-1.55 


2.92 


mW m^-^ 


0.81 


loglo(M) 


0.89 


0.11 


0.01 


-0.15 


-0.07 


0.04 




0.82 


^procl^s7n 


1.04 


0.00 


0.06 


-0.01 


0.07 


0.16 




0.74 


logio(Aic,tot) 


8.42 


-0.15 


-0.16 


0.22 


0.20 


0.01 




0.87 



Table 4: Coefficients for the best-fit linear function (Eq. |22| relating parameter values after 4.5 Gyr for parameters with correlation 
coefficients > 0.70 to a given set of initial conditions for Venus. Correlation coefficients quantifying the correspondence between the 
actual and predicted output parameters were calculated using normalized and mean subtracted input and output parameters. The 
average values of the input parameters are Tu(0) = 1700 K, Tcm(O) = 4000 K, log^oCw) = 19i logio(^^™) ~ 1' ^'^'^ Qo = 1-05 X lO"'^ 
W m~^. For the best-fit function, the input parameters are mean subtracted and normalized by 212 K, 408 K, 0.82, 0.82, and 4.30 X 
10~* W m""^, respectively. 
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Table 5: Coefficients for the best-fit linear function (Eq. \22\ relating parameter values after 4.5 Gyr for parameters with correlation 
coefficients > 0.70 to a given set of initial conditions for a 5M(^ super- Venus. Correlation coefficients quantifying the correspondence 
between the actual and predicted output parameters were calculated using normalized and mean subtracted input and output parameters. 
The average values of the input parameters are Tu(0) = 1701 K, Tcm(O) = 4351 K, log]^Q{r7o) = 19i logioC^^m) = li and Qo = 1-73 X 
10~^ W m~''. For the best-fit function, the input parameters are mean subtracted and normalized by 212 K, 408 K, 0.82, 0.82, and 7.07 
X 10~* W m""^, respectively. 



crust. From Eq. 12 volumetric melt production for a 
planet may scale as 



more than cores with increasing planetary mass. Next, 
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where the subscript © denotes values for an Earth-mass 
planet and Am stands for the mantle surface area. 

We can approximate S using the representative in- 
terior models of Valencia et al. ( 2006 1 , for which R ex 
P„ oc M^^'^^^Ta^^^dJocM^. First, consider 
the thickness of a melting region, dm = Zi - Zf. Since Zf 
= Pf/{pLg) is approximately constant for any planet, 
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where a roughly constant mantle to core thickness ratio 
is assumed, although planetary mantles grow slightly 
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The Rayleigh number for a massive planet scales as 
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Assuming that the first and second terms on the right 
hand side are roughly equal to unity, we have 



and thus 
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Figure 7: Sample histories for IMg) (solid lines), SMgj (dashed lines), and IOMq (dotted lines) super-Venus planets. Red, blue, and 
black curves, respectively, signify (a) crust, mantle potential, and core/mantle boundary temperatures, (b) surface, mantle, and core heat 
flows, (c) crust, depleted mantle lithosphere, and mantle lithosphere thicknesses, and (d) normalized mantle water content and fraction 
of processed source mantle. Default initial conditions are Qo = 1.0 xlO"^ W m-3 (scaled with pm), T„(0) = 1700 K, rjo = 10^^ Pa ■ s, 
Ar/™ = 100, and {dp/d<t>) = 120 kg m-3. The 1, 5, and lOM^ planets have Tc™,(0) = 4000, 4350, and 4900 K, respectively. Because 
crustal melting causes highly discontinuous surface and mantle heat fluxes, a moving average with a 75 Myr span was used for plotting 
purposes. 



Because he is usually much smaller than Rp 

A.rn ( Rp — \ / Rp 



R 



M 



(29) 

Finally, the rest of the scaling relations may simply be 
assumed as 



1 



and 



h"m 



M 
Mm 



(30) 

62 

(31) 

Hence, S « 0.240 and (/™//™,©) « (M/Afe)"''""- 
Because he ~ fm ^ At/^AnRp), where At is the du- 
ration of crust growth, an increase in melt produc- 
tivity with mass does not guarantee an increase in 
crustal thickness with mass. As planetary mass, and 
thus radius, increases, a larger volumetric melt produc- 
tion is required to produce a certain crustal thickness. 
Specifically, crustal thickness would only increase with 
mass for 6 > 0.524, assuming that At is roughly con- 



stant. Therefore, although melt productivity increases 
with planetary mass, this simple scaling analysis indi- 
cates that crustal thickness should decrease with scaling 



-0.284 



Panel (a) of Fig. 10 is a plot of model output present- 
day crustal thickness as a function of planetary mass 
for simulations of Mars, Venus, and seven super- Venus 
planets. While initial conditions strongly affect simu- 
lation results, the model outputs generally follow this 
simple scaling. Smaller planets can have thicker crust 
though they tend to be characterized by lower mantle 
temperatures. 

^.^.2. Mantle Processing 

The scaling of mantle processing with planetary mass 
follows easily from the above analysis. A simplified 
equation for the volume of processed mantle is 



fn 



V sa ^At 



(32) 



where At is a duration for crustal growth. 

Thus, the amount of processed mantle scales with 
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Figure 8: Summary of parameter values at the present for 672 simulations of the thermal evolution of a 5Mq super- Venus. Arrows are 
projections of the principal component basis vectors that emanate from a point representing the averaged simulation results, indicating 
axes that account for the vast majority of the data set's variance. The red arrow represents a larger percentage of cumulative variance 
(41%) than the green arrow (28%). Panels show (a) Moho temperature, (b) mantle potential temperature, (c) surface heat flux, (d) 
mantle heat flux, (e) fraction of mantle processed by melting, and (f ) total time for crust to grow from 10% to 95% of its present thickness 
as functions of crustal thickness. Because crustal melting causes highly discontinuous surface and mantle heat fluxes, the model output 
is the averaged values for the flnal 100 Myr of planetary evolution. 



planetary mass as 



^proc 



fr, 



M 
Mm 



(33) 



so^^Sk 0.240. 

The volume of a super- Venus planet scales as 



V 



R 



M 

Mm 



(34) 



-0.546 



SO C = 0.786. Therefore, (Vproc/V) cx {M/Mfs) 
Although the amount of processed mantle material in- 
creases with planetary mass, the fraction of processed 
mantle decreases with increasing planetary mass be- 
cause the mantle volume increases more rapidly than 
the amount of processed material. Panel (b) in Fig. 10 
confirms that the fraction of processed mantle does in- 
deed decrease with increasing planetary mass according 
to this scaling law, although initial conditions strongly 
affect the simulation results. 

4-5. Viscosity Contrasts During Stagnant- Lid Convec- 
tion 

The viscosity contrast across the lithosphere is 
tracked during each thermal evolution simulation, along 
with the critical viscosity contrast above which a planet 
is locked in the stagnant-lid regime. Figure 11 shows 



the output of 595 simulations for Mars, Venus, and two 
super- Venus planets for which Qg, T„(0), and /i were 



varied over a wide range. In particular, all permuta- 
tions of Qo = 0.5, 1.0, and 1.75 xlO"'^ W m^^ (scaled 
as usual with p,„) and r„(0) = 1400, 1700, and 2000 K 
were considered for a range of fi between 0.0 and 0.9. 

From this plot, several conclusions may be drawn. 
First, for values of the frictional coefficient associated 
with dry silicate rocks, fi ^ 0.7 to 0.8, plate tectonics is 
never favored. Second, increasing planetary mass does 
not substantially affect the likelihood of plate tectonics. 
Third, the effects of choosing different initial conditions 
are amplified for greater planetary mass. Finally, al- 
though choosing extreme initial conditions can change 
the viscosity contrast by orders of magnitude, the effect 
of the friction coefficient is far more important. 



4.6. Formation of an Eclogite Layer 

At depth, crustal rock may undergo a phase transi- 
tion to eclogite. To extend the simple analysis from 
earlier, we write the thickness of the crust in the eclog- 
ite stability field as = he - de, where de is the depth 
of the phase boundary. Likewise, we consider 6^ and 6a, 
the fractions of the crust in and above, respectively, the 
eclogite stability field. Because (5e+(5a = Se.e+^a,® = 1: 
we may write 



1 



Sa 



(35) 
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Figure 9; Summary of parameter values at the present for 675 simulations of the thermal evolution of a lOMg super- Venus. Arrows are 
projections of the principal component basis vectors that emanate from a point representing the averaged simulation results, indicating 
axes that account for the vast majority of the data set's variance. The red arrow represents a larger percentage of cumulative variance 
(42%) than the green arrow (30%). Panels show (a) Moho temperature, (b) mantle potential temperature, (c) surface heat flux, (d) 
mantle heat flux, (e) fraction of mantle processed by melting, and (f ) total time for crust to grow from 10% to 95% of its present thickness 
as functions of crustal thickness. Because crustal melting causes highly discontinuous surface and mantle heat fluxes, the model output 
is the averaged values for the flnal 100 Myr of planetary evolution. 



The fraction of the crust above the eclogite stabihty 
field may scale as 



del he 



5aA 



dp 



M 

Mm 



(36) 

If we assume that pressure increases hydrostatically 
with depth and that the critical pressure below which 
the phase transition occurs is a constant, then d^ ~ l/g. 
So, e = 0.284 - 0.503 = -0.219. Therefore, the fraction of 
crust in the eclogite stability field should increase with 
planetary mass as 



1 - 



M 



-0.219 



(37) 



In thermal evolution models, the heat conduction 
equation is numerically solved to calculate crustal tem- 
peratures. An approximate temperature profile can 
also be calculated using a steady-state approximation 



as (Turcotte and Schubert 2002) 



T(z) = T +—Z - 9Ahlz^ 



(38) 



where Qc, the volumetric crustal heat production, is 
calculated as 



Qdtp) = Qoe 



(39) 



where tp is 4.5 Gyr and Vc is the volume of the crust. 
The boundary condition T{hc) = Tc is used to calculate 
surface heat flux for specified Moho and surface temper- 
atures and magnitude of internal heat production. Fi- 
nally, Eq. [38]is used to calculate the temperature profile 
throughout the entire thickness of the crust. Represen- 
tative temperature profiles for a planet can be used to 
approximate the fraction of crust that lies within the 
eclogite stability field. 



Figure [12] shows representative temperature profiles 
for Venus and 5 and lOM® super- Venus planets, calcu- 
lated using representative crustal thicknesses, degrees 
of mantle processing, and Moho temperatures from the 
previous sensitivity analyses. A range of internal radio- 
genic heating was also considered. The stability field of 



eclogite is taken from Philpotts and Ague (2009) and 



is drawn assuming a hydrostatic pressure increase with 
depth. Panel (d) in Fig. 12 summarizes the effects of 



initial conditions on the fraction of crust in the eclog- 
ite stability region after 4.5 Gyr and shows the scaling 
from Eq. |37j For Mars, Venus, and seven super- Venus 
planets, 54 thermal evolution simulations were run to 
study all permutations of the initial conditions r„ (0) ~ 
1700 and 2000 K and Venus-equivalent Qo = 0.5, 1.0, 
and 1.75 xl0~^ W m^"^. As predicted, the fraction of 
eclogite crust increases with planetary mass. 
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Figure 12: Crustal temperature profiles for (a) Venus and (b) and (c) IOMq super- Venus planets, calculated assuming representative 

crustal thicknesses, degrees of mantle processing, and Moho temperatures. The green shaded area is the approximate stability field for 
eclogite, drawn using the phase diagram from |Philpotts and A gue (2009J and assuming a hydrostatic pressure increase with depth. Black 
solid and red dashed lines represent Venus-equivalent Qo = 1-0 x 10~^ and 1.75 X 10~^ W m""^, respectively. Panel (d) shows the 
fraction of crust in the eclogite phase for Mars, Venus, and seven super- Venus planets. Circles and triangles represent r„(0) = 1700 and 
2000 K, respectively. Blue, black, and red symbols represent Venus-equivalent Qo = 5-0 x 10~^, 1.0 X 10~^, and 1.75 X 10~^ W m"'', 
respectively. 



5. Discussion 

5.1. Pressure Effects on Mantle Rheology 

The rheological behavior of the mantles of large 
rocky planets is difficult to predict. While the 
core/mantle boundary pressure for Earth is ^135 GPa 
(e.g., Dziewonski and Anderson 1981 1, pressures within 
the silicate mantles of large rocky planets likely exceed 
1 TPa 



(.e.g., 



Valencia et al. 2006). Above the tran- 



sition to post-perovskite at ^120 GPa, Earth's man- 
tle is primarily made of MgSiOa-perovskite and MgO 
(e.g., Murakami et al. 2004). In contrast, much of 



the mantles of super- Venus planets will be dominated 
by post-perovskite and perhaps, above 1 TPa, a mix- 
ture of MgO and SiOa (lUmemoto et all |2006p. Un- 



fortunately, we lack experimental measurements of the 
properties of planetary materials under these extreme 
conditions. Until such data are available, conjectures 
about the rheology of the silicate mantles of large rocky 
planets will remain controversial. Thermal evolution 



simulations are very sensitive to assumed rheological 
behaviors, so investigating the implications of various 
possible assumptions is essential. 

The viscosities of most planetary materials increase 
with pressure when examined at relatively low pressures 



(e.g., jKaxato 2008). Simple extrapolation of this trend 
predicts extreme increases in viscosity within massive 
terrestrial planets. Extensions of known perovskite rhe- 
ology, for instance, imply an increase of >15 orders of 
magnitude as pressure increases to 1 TPa in an adia- 



batic mantle (Stamenkovic et al. , 2011). Specifically, 
a viscosity profile may be calculated as ( [Stamenkovic| 



et al. 2012) 



r]{P,T) = ?7oexp 







1-] 








[t 


T* J 


4 





, (40) 



where rjQ is a reference viscosity at the reference tem- 
perature T* — 1600 K and V* is an activation volume. 
Figure [13] shows calculated viscosity profiles for rjo = 
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Figure 11: Summary of 595 simulations of the evolution of Mars, 
Venus, and two super- Venus planets, showing the minimum ra- 
tio of actual viscosity contrast to critical viscosity contrast and 
thus the likelihood of plate tectonics being favored at some point 
during 4.5 Gyr of planetary evolution. Each planet was evolved 
from six different sets of initial conditions (three values for both 
radiogenic heating and mantle potential temperature) for many 
different values of fi, the effective friction coefficient. Points plot- 
ted above the indicated line represent simulations for which the 
actual viscosity contrast never dipped below the critical value for 
a transition to plate tectonics. Below the indicated line, which 
occurs only for fi < 0.3, plate tectonics may have been favored at 
some point. For dry silicate rocks, fi ~ 0.7 to 0.8. 



Figure 10: Summary of 54 simulations of the evolution of Mars, 
Venus, and seven super- Venus planets, showing the correspon- 
dence between simulation results and simple scaling laws for the 
effects of planetary mass on (a) crustal thickness and (b) mantle 
processing. Circles and triangles represent Tu(0) = 1700 and 2000 
K, respectively. Blue, black, and red symbols represent Venus- 
equivalent Qo = 5.0 X 10-8, 1.0 X lO"'^, and 1.75 x 10"'^ W 
m~^, respectively. Dashed black lines show the scaling relations 
(a) he oc {Af/Me)-''-2«4 and (b) (Vproc/V) oc {M / M(q)-°-^^<^ , 
with each curve fixed to intersect the average output from the 
simulations for the 2M0 super- Venus planet. 



10^1 Pa s and V* = 2.5, 1.7, and 0.0 cm^ mori for 
the convecting, adiabatic mantle within a IQM^ super- 



60 v'' = 2.5cm^ mo!"'' 



Venus planet, following iStamenkovic et al. (2012). Our 



parameterized formulation for stagnant-lid convection 
is based on numerical modeling with the incompress- 
ible fluid approximation using temperature-dependent 



but pressure-independent viscosity (e.g., Solomatov and 



Moresi 


2000 


Korenaga 


2009 



Therefore, we are as- 
suming pressure-independent constant potential viscos- 
ity, where the effect of temperature increase along an 
adiabatic gradient exactly balances the effect of pres- 
sure on viscosity, which requires a non-zero, positive 
activation volume. With the linear temperature gra- 
dient in Fig. |13[ the assumption of constant potential 
viscosity corresponds to V* — 0.22 cm'^ mol~"'^. 

Strongly pressure-dependent viscosity can cause dra- 
matically different behavior to emerge from parameter- 
ized convection models, including sluggish lower mantle 
convection and even the formation of a conductive lid 
above the core/mantle boundary (IStamenkovic et al 
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Figure 13: Internal temperature and viscosity profiles for a lOMgj 
super- Venus planet, calculated as in [Stamenkovic et al.| | |2012 l. 
The black curve shows internal temperature as a function of depth 
in the convecting mantle. Green, red, and blue dashed lines rep- 
resent viscosity profiles calculated using Eg. |40| for V* = 0, 0.22, 
1.7, and 2.5 cm^ mol~^, respectively. 



lower mantle of large rocky planets, melt production 
would be significantly decreased and the likelihood of 
plate tectonics might decrease along with convective 
vigor as planetary mass increased. The thermal conduc- 
tivity and expansivity of the mantle are also predicted 
to increase and decrease, respectively, with depth be- 
cause of increasing pressure, but the effects of these 
changes on mantle dynamics are dwarfed by the pu- 



2012). If convection were effectively suppressed in the tative increase in mantle viscosity (Stamenkovic et al 
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2011 2012) 



Although an increase in viscosity under greater pres- 
sures seems intuitive, the straightforward apphcation of 
hmited, low-pressure experimental data may not accu- 
rately describe the rheology of massive terrestrial plan- 
ets. In fact, four mechanisms, including a transition 
from vacancy to interstitial diffusion mechanisms, may 
cause a viscosity decrease with depth above a pres- 
sure of ^0.1 TPa (Karato 2011), as post-perovskite 



On Earth, the phase transition from (metamor- 
phosed) basalt to eclogite primarily occurs in sub- 
duction zones. Hydration may thus be important 
to allowing this phase transition to occur relatively 



rapidly (e.g., Ahrens and Schubert 1975 ), although this 
type of kinetic calculation strongly depends on diffu- 



and additional high-pressure mineral phases dominate 
mantle rheology. According to this study, viscosities 
in the deep interiors of super- Venus planets may be 
less than the viscosity of Earth's lower mantle, poten- 
tially by as much as 2-3 orders of magnitude. More- 
over, the depth-dependence of viscosity within Earth's 
perovskite-dominated mantle is still debated because 
Earth's viscosity profile has not been well-constrained. 

Despite basic consensus that Earth's lower mantle is 
probably more viscous than the upper mantle, the mag- 
nitude of the viscosity contrast remains controversial. 
Early studies of Earth's topography and geoid suggested 
a ~300-fold increase in viscosity between Earth's up- 
per and lower mantle (Hager and Richards 1989). On 



sion data that are not well-constrained (e.g., Namiki 
and Solomon 1993). Eclogite is also formed during 
continent-continent collisions such as the Eurasian and 
Indian plate collisions ( Bucher and Prey 2002 1 . Fur- 



thermore, the high density of eclogite is theorized to 
have caused delamitation, foundering, and recycling 
of relatively thick oceanic lithosphere on Earth dur- 
ing the Archaean (Vlaar et al. 1994). Finally, eclog- 



ite may be produced in large mountain ranges through 



magmatic differentiation (Ducea 2002) and pressure- 
induced phase transition ( Sobolev and Babeyko^ 2005| 
in thick continental crust. Evidence for the strong influ- 
ence of recent eclogite production and foundering on the 
topography of the central Andes Mountains has been 
gathered through geodynamics, petrology, and seismol- 
ogy (e.g., Kay and Abbruzzi', 1996 Beck and Zandt 



the other hand, analyses of post-glacial rebound predict 
an order of magnitude less of viscosity increase (e.g., 
Kaufmann and Lambeck 2002) and gravity data are 



as synthesized in a numerical study (Pelletier et al 



2002[ [Sobolev and Babeyko , 2005[ [Schurr et al.[|2006[). 



consistent, albeit loosely, with uniform or only slightly 



depth-dependent mantle viscosity (Soldati et al. 2009). 
A joint inversion of these data sets predicts that Earth's 
internal viscosity increases by ~2 orders of magnitude 



throughout the mantle (Mitrovica and Forte 2004), but 



2010). Because massive terrestrial planets have rela- 
tively high surface gravity, the phase transition to eclog- 
ite will occur at a comparatively shallow depth, making 
eclogite the stable mineral phase for a large fraction 
of the crust. The formation of a thick eclogite layer 
then could cause lithosphere foundering or intermittent 
plate tectonics, as has been proposed in episodic sub- 



it is noted that such inversions are known to suffer duction mechanisms for Venus (Turcotte, 1993 Fowler 



King 


1995 


Kido and 


and O'Brien 


1996 



Cadek 1997). In any case, considering the tremendous 
uncertainty surrounding viscosity profiles within Earth 
and putative super- Venus planets, it remains legitimate 
to investigate how large rocky planets might evolve in 
the limiting case of pressure-independent potential vis- 
cosity. 

5.2. Escaping the Stagnant-Lid Regime 

Terrestrial planet evolution strongly depends on the 
regime of mantle convection. Assuming that brittle fail- 
ure limits the strength of the lithosphere, our simula- 
tions indicate that the effects of lithosphere hydration 
dominate the effects of planetary mass on yield and 
convective stresses. That is, the increase in convec- 
tive vigor with planetary mass only makes plate tecton- 
ics marginally more likely. Modeling results for super- 
Venus planets, however, suggest two additional mecha- 
nisms for escaping the stagnant-lid regime. First, mas- 
sive terrestrial planets in the stagnant-lid regime feature 
crustal temperature profiles that enter the stability field 
of eclogite after crust grows beyond a critical thickness. 
If a sufficiently large fraction of the total crustal thick- 
ness is composed of eclogite, the entire crust could be 
gravitationally unstable and susceptible to foundering 
because eclogite is intrinsically denser than mantle peri- 
dotite. 



Representative temperature profiles for massive ter- 
restrial planets pass through the eclogite stability field 
for plausible initial conditions. In fact, radiogenic heat- 
ing and thus crustal temperatures should be greater 
than calculated with Eq. [38j which would increase 
the speed of the phase transition to eclogite, because 
crustal construction mostly occurs early in planetary 
history. For super- Venus planets with masses greater 
than ~4M0, eclogite may be the stable phase for the 
majority of the crust unless the initial mantle potential 
temperature or magnitude of internal heating is very 
low. So, crust material may undergo a phase tran- 
sition to eclogite at relatively shallow depths as the 
crust grows during thermal evolution in the stagnant-lid 
regime, forming a thick eclogite layer that could subse- 
quently founder. The buoyant stress from the presence 
of eclogite scales as Apgh^, where Ap ~ 100 kg m~^ is 
the difference between the densities of eclogite and the 
mantle. On the other hand, the lithospheric strength 
scales as fj,pLgde^ which can only be overcome when the 
depth scale of the eclogitic layer becomes large enough 
(at least locally, for example, by foundering). As long 
as crustal production continues on large rocky planets, 
eclogite formation and foundering could occur period- 
ically, possibly yielding a regime of mantle convection 
resembling intermittent plate tectonics. Although we 
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suggest that this process is plausible, pursuing its dy- 
namics in detail is left for future studies. 

High surface and crustal temperatures may also cause 
periodic transitions from the stagnant-lid regime to a 
form of mobile-lid convection. In this work, massive ter- 
restrial planets in the stagnant-lid regime with surface 
temperatures held constant at 300 K tend to have very 
hot crusts. If high surface temperatures exist alongside 
high crustal temperatures, a transition from stagnant- 



lid convection to a mobile- lid regime can occur (Reese 
et al. 1999| ) . Feedback between a changing mantle con- 
vection regime and a periodic atmospheric greenhouse 
effect driven by varying amounts of volcanism, for in- 
stance, may be very important to the evolution of Venus 
(Noack et al., 2011). As surface temperature depends 



on atmospheric mass and the composition and lumi- 
nosity of the central star, however, this possibility of 
escaping the stagnant-lid regime may not be as robust 
as the first mechanism based on the formation of self- 
destabilizing crust 

5.3. Limitations of Parameterized Models 

Any parameterized model suffers shortcomings. 
Steady-state evolution is assumed, for instance, which 
poorly captures transient events that occur early in 



planetary evolution such as large impacts ( Agnor et al 



19991 and the crystallization of a magma ocean (e.g., 
Solomatov and Stevenson 1993). Fundamental as- 



sumptions such an adiabatic temperature gradient in 
the mantle and pressure-independent potential viscos- 
ity are controversial, and different approaches such as 
mixing length theory (I Wagner et al. 2011 ) may be nec- 



essary to calculate the thermal structure of planetary 
interiors if they are not valid. However, our simpli- 
fied simulations only aim to illuminate the first-order, 
relative effects of planetary mass on terrestrial planet 
evolution. Recreating the thermal history of a particu- 
lar planet would require the introduction of many addi- 
tional complications. One-dimensional models only re- 
turn globally averaged values for important quantities, 
for instance, but mantle plumes, which may upwell from 
the core/mantle boundary when the core heat flux is 
positive, are likely important to local magmatism and 
surface features on terrestrial planets hke Mars (e.g., 
Weizman et al. 2001 1 and Venus (e.g., Smrekar and 



Sotin 2012). Furthermore, applying a parameterized 



approach to Venus, where the precise quantity of mag- 
matism is an key output, requires more computationally 
intensive simulations to benchmark the relevant scaling 
laws. 



yield first-order, relative conclusions about the evolu- 
tion of generic terrestrial planets in the stagnant-lid 
regime. Principal component analysis of simulation re- 
sults conducted with a wide range of initial conditions 
captures the relationships between the large number of 
parameters that describe the interior of a planet. De- 
pending on initial conditions, these planets may have 
evolved along a variety of paths, featuring different 
crustal thicknesses and temperatures, interior tempera- 
tures, and degrees of mantle processing. To produce 
specific histories consistent with spacecraft data ob- 
tained from Mars and Venus, complications must be 
added to these simple models. 

Properties of massive terrestrial exoplanets are poorly 
constrained, so questions about the effects of planetary 
mass on the likelihood of plate tectonics and other im- 
portant planetary parameters await definitive answers. 
In this study, we explored what might happen if in- 
ternal viscosity is not strongly-pressure dependent, the 
alternative to which has been explored previously us- 
ing parameterized models. Although convective vigor 
increases with planetary mass, the likelihood of plate 
tectonics is only marginally improved. Simple scaling 
analyses indicate that mantle melt productivity should 
increase with planetary mass. Because the increase in 
mantle processing is slow, however, crustal thickness 
and the relative fraction of processed mantle actually 
decrease with increasing planetary mass, as thermal 
evolution simulations confirm. Surface gravity increases 
with planetary mass, so pressure in the crust of mas- 
sive terrestrial planets increases relatively rapidly with 
depth. Plausible temperature profiles favor a phase 
transition to gravitationally unstable eclogite during 
normal crustal formation, whereas the basalt to eclog- 
ite transformation rarely occurs aside from subduction 
on Earth. Therefore, thick eclogite layers, along with 
mobile, hot crustal material, may be important to the 
evolution of massive terrestrial planets. 
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Table A.l: Principal component basis matrix for Venus (V) and a IOM0 super- Venus planet (S) for the model output after 4.5 Gyr of 
planetary evolution. Two eigenvectors account for over 65% of the variance in the normalized and mean subtra<;ted simulation results. 
The fractions of the cumulative variances for which each principal component accounts, calculated by dividing the principal component 
eigenvalue by the sum of the eigenvalues for all principal components, are in the bottom row. Output parameters were mean subtracted 
and normalized using the listed average and standard deviation values. 



Appendix A. Statistical Analysis of Simulation 
Results 

Parameterized evolution models involve quite a few 
model parameters. It is important to understand how 
simulation results depend on a particular choice of 
model parameters by testing a variety of situations, but 
at the same time, it becomes difficult to grasp the in- 
flated amount of numerical data. Principal component 
analysis (PCA) can be used to assess the effective di- 
mensionality of a given data space. Our intention here 
is to use PCA to extract major features and trends from 
a large number of simulation results. Each sensitivity 
analysis consists of n simulations with m output pa- 
rameters, comprising a data set D™. Some parameters, 
such as Af]^ and u, exhibit orders of magnitudes of vari- 
ation, and we consider their logarithms because PCA is 
designed for linear data sets. We normalize the data set 
as 



pm _ 



(A.1) 



where /x™ is the average value of the m-th output pa- 
rameter. 



1 



(A.2) 



which account for a progressively decreasing percentage 
of data variance. Principal components accounting for 
at least 65% (an arbitrary threshold) of the total vari- 
ance are selected for examination to reveal important 
aspects of simulation results. 



and a™ is the standard deviation of the m-th output 
parameter. 



1 " 
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(A.3) 



Because the normalized data have zero mean, the 
covariance matrix Cp = P can be decomposed as 
Cp = ■ diag[Ai . . . Am] • A, where A^, the eigenvalues, 
are ordered so that Ai > A2 > . . . > Xm- The cor- 
responding eigenvectors are the principal components. 
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